// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#ifndef ffFUNCTION_H
#define ffFUNCTION_H

#include "type.h"

//! Some neccessary function.
namespace function
{
	const double pi = 3.14159265358979323846;            //! pi
	//! Convert deg to radian.
	inline double deg2rad(const double & deg );
	//! Convert radian to deg.
	inline double rad2deg(const double & rad ) ;

	inline atom_index       ait2ai(array_index_type i);
	inline array_index_type ai2ait(atom_index i);

	//! Sort and swap. As result: a < b
	template<typename T>
	inline void sort(T &a,T &b) { if(b < a) std::swap(a,b); }

	//! Sort and swap three variables. As result: a < b < c
	template<typename T>
	inline void sort(T &a,T &b,T &c)
	{ if(b < a) std::swap(a,b); if(c < a) std::swap(a,c); if(c < b) std::swap(b,c);	}
	
	//! Calculate distance of two vector.
	inline double dist(const vector_type & a,const vector_type & b);
	//! Calculate length of a vector.
	inline double length(const vector_type & a);
	//! x^2
	inline double sqr(const double & x);
	//! Calculate cos(a) of two vector.
	inline double cosa(const vector_type & a,const vector_type & b);

	inline void cos_angle_deriv(const vector_type & a,const vector_type & b,vector_type & da,vector_type & db);
	inline void cos_dihedral_deriv(const vector_type & a,const vector_type & b,const vector_type & c,
								  vector_type & da,vector_type & db,vector_type & dc); 
	inline void cos_inversion_deriv(const vector_type & a,const vector_type & b,const vector_type & c,
								  vector_type & da,vector_type & db,vector_type & dc);
}


//! Convert deg to radian.
inline double function::deg2rad(const double & deg ) { return deg*pi/180.00; }
//! Convert radian to deg.
inline double function::rad2deg(const double & rad ) { return rad*180.00/pi; }

inline atom_index       function::ait2ai(array_index_type i)  { return i+1; }
inline array_index_type function::ai2ait(atom_index i)        { return i-1; }

//! Calculate distance of two vector.
inline double function::dist(const vector_type & a,const vector_type & b) { return (a-b).norm(); }
//! Calculate length of a vector.
inline double function::length(const vector_type & a) { return a.norm(); }
//! x^2
inline double function::sqr(const double & x) { return x*x; }
//! Calculate cos(a) of two vector.
inline double function::cosa(const vector_type & a,const vector_type & b) { return a.dot(b)/a.norm()/b.norm(); }

inline void function::cos_angle_deriv(const vector_type & a,const vector_type & b,vector_type & da,vector_type & db) 
{ 
	real_type la = a.norm();
	real_type lb = b.norm();
	real_type alpha = 1.0/(la*lb);
	real_type beta = lb/la;

	da = alpha*b - alpha*alpha*beta*a;
	db = alpha*a - alpha*alpha/beta*b;
}

inline void function::cos_dihedral_deriv(const vector_type & a,const vector_type & b,const vector_type & c,
								  vector_type & da,vector_type & db,vector_type & dc) 
{ 
	vector_type axb = a.cross(b);
	vector_type bxc = b.cross(c);
	vector_type bxbxc = b.cross(bxc);
	vector_type bxaxb = b.cross(axb);
	vector_type cxaxb = c.cross(axb);
	vector_type axbxc = a.cross(bxc);
	vector_type axaxb = a.cross(axb);
	vector_type cxbxc = c.cross(bxc);

	real_type laxb = axb.norm();
	real_type lbxc = bxc.norm();

	real_type alpha = 1.0/(laxb*lbxc);
	real_type beta = lbxc/laxb;

	da = alpha*bxbxc - alpha*alpha*beta*bxaxb;
	db = alpha*(cxaxb-axbxc) - alpha*alpha*(cxbxc/beta-beta*axaxb);
	dc = -alpha*bxaxb + alpha*alpha/beta*bxbxc;
}

inline void function::cos_inversion_deriv(const vector_type & a,const vector_type & b,const vector_type & c,
								  vector_type & da,vector_type & db,vector_type & dc) 
{ 
	vector_type axb = a.cross(b);
	vector_type bxc = b.cross(c);
	vector_type axc = a.cross(c);
	

	real_type laxb = axb.norm();
	real_type lc = c.norm();

	real_type alpha = 1.0/(laxb*lc);
	real_type beta = lc/laxb;

	da = alpha*bxc - alpha*alpha*beta*b.cross(axb);
	db = -alpha*axc + alpha*alpha*beta*a.cross(axb);
	dc = alpha*axb + alpha*alpha/beta*c;
}




#endif // ffFUNCTION_H

